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Abstract 

A  global  barotropic  model  of  the  atmosphere  is  presented  governed  by  the  shallow 
water  equations  and  discretized  by  a  Runge-Kutta  discontinuous  Galerkin  method 
on  an  unstructured  triangular  grid.  The  shallow  water  equations  on  the  sphere,  a 
two-dimensional  surface  in  M3,  are  locally  represented  in  terms  of  spherical  triangu¬ 
lar  coordinates,  the  appropriate  local  coordinate  mappings  on  triangles.  On  every 
triangular  grid  element,  this  leads  to  a  two-dimensional  representation  of  tangential 
momentum  and  therefore  only  two  discrete  momentum  equations. 

The  discontinuous  Galerkin  method  consists  of  an  integral  formulation  using  a 
Rusanov  numerical  flux.  A  strong  stability-preserving  third  order  Runge-Kutta 
method  is  applied  for  the  time  discretization.  The  polynomial  space  of  order  k  on 
each  curved  triangle  of  the  grid  is  characterized  by  a  Lagrange  basis  and  requires 
high-order  quadature  rules  for  the  integration  over  elements  and  element  faces.  For 
the  presented  method  no  mass  matrix  inversion  is  necessary,  exept  in  a  preprocessing 
step. 

The  validation  of  the  atmospheric  model  has  been  done  considering  steady-state 
and  unsteady  analytical  solutions  of  the  nonlinear  shallow  water  equations.  Exper¬ 
imental  convergence  was  observed  and  the  order  of  convergence  k  +  1  was  achieved. 
Furthermore,  the  article  presents  a  numerical  experiment,  for  which  the  third  order 
time-integration  method  limits  the  model  error.  Thus,  the  time  step  At  is  restricted 
by  both,  the  CFL-condition  and  accuracy  demands.  As  a  second  step  of  valida¬ 
tion,  the  model  could  reproduce  a  known  barotropic  instability  caused  by  a  small 
initial  perturbation  of  a  geostrophic  balanced  jet  stream.  Conservation  of  mass  was 
shown  up  to  machine  precision  and  energy  conservation  converges  with  decreasing 
grid  resolution  and  increasing  polynomial  order  k. 
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1  Introduction 


Modeling  atmospheric  flows  for  climate  simulations  as  well  as  for  weather  prediction 
is  a  complex  problem,  clue  to  the  nonlinear  structure  of  the  dynamical  and  physi¬ 
cal  phenomena  on  widely  varying  spatial  and  temporal  scales  and  their  multi-scale 
interaction  processes.  Depending  on  the  complexity  of  an  atmospheric  model  the 
governing  equations  are  the  fundamental  atmospheric  conservation  laws  for  mass, 
momentum  and  energy  or  appropriate  simplifications  of  them.  If  the  regarded  equa¬ 
tion  set  is  a  hyperbolic  system,  energetic  shocks  can  develop  theoretically.  Although 
this  is  usually  not  the  case  in  atmospheric  models,  the  discretization  should  repre¬ 
sent  regions  of  scale  collapse  and  breaking  waves  generating  discontinuities  in  the 
velocity  field;  the  discrete  conservation  properties  of  the  DG  method  are  appropriate 
for  this  task. 

The  shallow  water  equations  (SWE),  valid  for  a  homogeneous  atmosphere  with 
small  vertical  velocities  and  horizontal  velocities  independent  in  the  vertical  direc¬ 
tion,  constitute  a  hyperbolic  system  of  conservation  laws.  It  is  one  of  the  simplest 
nonlinear  hyperbolic  systems,  covering  important  planetary  atmospheric  features, 
like  the  Rossby  wave  formation. 

For  the  spherical  SWE  the  spatial  domain  is  the  sphere  S,  a  two-dimensional 
surface  in  M3.  In  a  regional  or  mesoscale  SWE  model  the  momentum  is  a  two- 
dimensional  vector.  In  contrast,  the  Cartesian  formulation  of  the  spherical  case  in 
[8]  represents  the  tangential  momentum  of  the  flow  as  a  three-dimensional  vector 
and  includes  a  Lagrangian  multiplier  to  constrain  the  momentum  to  be  tangential. 
Applying  this  form  to  a  numerical  model  usually  leads  to  three  momentum  equa¬ 
tions  and  requires  a  correction  step  to  realize  the  constraint  discretly,  see  e.  g.  [13]. 
Models  in  standard  spherical  coordinates  realize  a  two-dimensional  momentum  rep¬ 
resentation  but  have  to  pay  additional  attention  to  phenomena  near  the  poles  due 
to  singularities  of  the  coordinate  mapping,  see  e.g.  [18]. 

The  idea  to  avoid  the  drawbacks  of  the  Cartesian  and  the  spherical  coordinates 
formulation  is  to  represent  the  spherical  SWE  in  terms  of  local  coordinate  trans¬ 
formations.  On  a  cubed-sphere  grid,  a  spherical  quadrilateral  grid,  [20]  and  [21] 
achieved  a  two-dimensional  momentum  representation  avoiding  any  pole  problem. 
The  presented  model  provides  these  properties  on  spherical  unstructured  triangular 
grids  using  spherical  triangular  coordinates.  Using  coordinate  independent  differ¬ 
ential  operators  on  S  known  from  differential  geometry,  the  spherical  SWE  are 
considered  in  flux  form  on  the  surface  S.  Spherical  triangular  coordinates,  that  is 
appropriate  local  coordinate  mappings  7 e  on  a  curved  triangle  E,  yield  the  two- 
dimensional  representation  of  tangential  momentum  vectors.  Numerous  numerical 
methods  have  been  proposed  for  future  global  atmospheric  models  including  finite 
volumes  [18],  [21],  spectral  elements  [25],  [10],  and  DG  methods.  We  have  selected 
the  DG  method  for  our  model  because  it  allows  us  to  achieve  high-order  accuracy 
as  in  spectral  elements  while  conserving  all  quantities  both  locally  and  globally  as 
in  finite  volumes.  Furthermore,  our  use  of  unstructured  triangular  grids  allows  for 
much  flexibility  in  future  work  on  adaptivity. 

DG  methods  have  been  already  successfully  applied  to  the  spherical  SWE  in  [12], 
[20]  and  [11].  DG  methods  in  combination  with  Runge-Kutta  (RK)  time  discretiza- 
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tions,  see  the  review  in  [5],  can  be  characterized  as  a  high-order  generalization  of 
finite  volume  methods.  This  method  is  robust,  high-order  accurate,  locally  conser¬ 
vative  with  high  potential  for  parallelization,  see  [11]. 

The  organization  of  this  article  is  as  follows.  In  Sec.  2  the  governing  spherical 
SWE,  are  given  using  surface  differential  operators.  Section  3  describes  the  numer¬ 
ical  discretization  by  a  RK-DG  method  applying  spherical  triangular  coordinates. 
In  Sec.  4  the  barotropic  model  of  the  atmosphere  based  on  the  discretization  is  val¬ 
idated  in  terms  of  steady-state  and  unsteady  analytical  solutions  and  a  barotropic 
instability  generated  by  a  small  initial  perturbation. 

2  Spherical  Shallow  Water  Equations 

The  spherical  shallow  water  equations  (SWE)  are  a  system  of  conservation  laws  for 
the  geopotential  layer  depth  (mass)  and  the  flow  momentum.  Because  the  integra¬ 
tion  domain  of  the  SWE  is  the  sphere,  a  two-dimensional  surface  in  M3,  the  system 
can  be  formulated  in  the  surrounding  Cartesian  space  M3,  see  [8].  Cote’s  formula¬ 
tion  is  equivalent  to  a  conservative  form  of  the  SWE  on  the  surface  S,  which  is  the 
formulation  used  further  below. 

Let  us  consider  the  sphere  S  =  {x  G  M  |  |x|  =  a}  with  the  Earth’s  radius  a  = 
6.371  x  106m,  the  geopotential  layer  depth  <f>  :  S  x  M+  — >  M,  the  tangential  velocity 
field  u  :  SxR+  — ■>  M3  with  u(x,  t )  G  TX(S )  and  the  conserved  variable  q  =  (<£>,  <l>uT)T. 
Then,  the  SWE  in  conservative  form  on  the  surface  S  are 

dtq  +  divs  f(q)  =  F(x,  q)  in  S  x  M+.  (1) 

Thereby,  the  flux  function  and  the  right  hand  side  are 

/(<?)  =  ’  /*(?)  =  fufa)  =$u®u  +  y Ids, 

F(x,  q)  =  iFu^x  q)  J  >:  Fu(x,  q)  =  ~fc$k  x  u  -  -  —k 

with  the  Earth’s  angular  velocity  f!  =  7.292  x  10~5s_1,  the  space  dependent  Coriolis 
parameter  fc(x)  =  212 ^j3-,  the  normal  unit  vector  k(x)  =  A  outward  on  S  and  the 
identity  mapping  M3  in  M3.  See  the  appendix  for  the  definition  of  the  differential 
operators  on  S.  Three  prognostic  equations  for  the  momentum  <f>u  appear  in  (1), 
whereas  the  momentum  is  forced  to  be  tangential  on  S  by  the  Lagrangian  multi¬ 
plier  —  ^F-k.  Due  to  this  forcing  term  no  global  conservation  of  momentum  can  be 
expected,  which  is  in  contrast  to  the  two-dimensional  shallow  water  equations.  As 
a  consequence  of  (1)  mass  is  locally  globally  conserved  while  energy  is  only  globally 
conserved. 


3  Discontinuous  Galerkin  Method 

The  Discontinuous  Galerkin  (DG)  method  is  applied  to  the  conservative  form  (1)  of 
the  SWE  on  the  surface  S.  Based  on  a  given  triangulation  on  every  curved  triangle 
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(element)  E  spherical  triangular  coordinates,  which  are  local  coordinate  mappings 
7 e  on  each  E,  are  introduced.  Using  7 E  polynomial  spaces  of  high-order  are  defined 
on  each  curved  element.  The  polynomial  representation  on  each  grid  element  uses 
high-order  Lagrange  polynomials  with  respect  to  Fekete  points.  This  approach  leads 
to  the  local  representation  of  the  tangential  momentum  fields  by  two  components 
only.  An  integral  form  of  (1)  leads  to  the  space-discrete  DG  method  including  a 
Rusanov  numerical  flux.  For  this  method  high-order  quadrature  rules  are  applied 
and,  except  in  a  preprocessing  step,  no  mass  matrix  inversion  has  to  be  evaluated. 
Finally  the  semi-discrete  problem  is  solved  by  a  strong  stability-preserving  explicit 
Runge-Kutta  (RK)  method.  The  fully  discrete  RK-DG  method  avoids  any  kind  of 
explicit  smoothing  such  as  diffusion  or  filter  operators. 


3.1  Spherical  Triangular  Coordinates 

Let  E  C  S  be  a  relative  open  spherical  triangle  bounded  by  great  circles  and  defined 
by  its  vertices  Xq,Xi,x2  €  S.  Then  we  define  for  E  the  local  coordinate  mapping 
7 e,  or  the  spherical  triangular  coordinates,  by 

je:D^E,  7E(9)=«g|jj.  (2) 

Here,  D  =  {y  G  M2  |  0  <  yi,  y2 ;  y\  +  y2  <  1}  is  a  two  dimensional  reference  triangle, 
and  xp(y)  =  x0+yi(xi—x0)+y2(x2—x0)  an  auxiliary  planar  mapping.  Following  the 
notation  in  the  appendix,  a  basis  of  TX(S),  Gram’s  determinant  and  the  Christoffel 
symbols  are  given  with  i,j,k  =  1,  2 


bi  = 


Xr, 


Xi 


X0) 


[Xi 


Xq) 


Xp  Xp 


Xr 


Xr, 


X. 


9  =  ip  1  X  b2)  ,  T)k  =  — |  •  {{Xj  -  x0)5lk  +  (xk  -  X0)6j). 


X, 


3.2  Discrete  Function  Space 

Let  T  =  {E  C  S\E  spherical  triangle,  E  open  in  S}  be  a  finite  conformal  trian¬ 
gulation  of  the  sphere  -  that  is,  a  triangular  grid  without  hanging  nodes.  For  the 
presented  model,  T  is  constructed  by  the  grid  generator  AMATOS,  see  [1].  Depen¬ 
dent  on  the  grid  level  /,  an  icosahedral  coarse  grid  Tq  is  refined  in  /  steps,  in  which 
every  triangle  of  T0  is  divided  by  bisection.  This  leads  to  an  unstructured  spherical 
triangulation  T  with  nearly  uniform  grid  resolution,  see  figure  1. 

The  polynomial  space  of  the  polynomials  of  degree  at  most  k  >  0  on  every  element 
EeT  is  defined  by 

P\E )  =  {<p:E^R\<polEe  Pk{D )}, 

where  7 E  is  the  spherical  triangular  coordinate  mapping  (2).  Thus,  the  coordinate 
mapping  7^  defines  as  well  the  curved  geometry  of  E  as  the  polynomial  space  on 
E.  This  technique  for  curved  elements  is  similar  to  the  definition  of  isoparametric 
finite  elements,  see  [3],  but  with  an  analytically  non-polynomial  mapping  ryE. 
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Figure  1:  Section  3.2,  Uniform  grid,  Grid  resolution  2058km  (left),  1041km  (middle), 
522km  (right). 


Every  polynomial  p  G  Pk(E)  is  represented  by  a  multivariate  Lagrange  basis 
(<Pi)i=i,..,Nk  of  Pk(E),  with  Nk  =  (/i+14A-'+2) ;  associated  with  the  Lagrange  points 
(xi)i=i,..,Nk)  that  is  for  x  G  E 


Nk 

P(x)  =  ^2<Pi(x)p(xi)-  (3) 

i— 1 

Thereby,  (ay)j=i,..,Ay  as  well  as  (</?j)*=i, ..,Nk  are  defined  by  tpi  =  &  o  y"1  and  xt  = 
7 E(xi).  (<Pi)i=i,..,Nk  is  the  multivariate  Lagrange  basis  of  Pk(D )  associated  with  the 
Fekete  points  (xi)i=i,..,Nk  in  D  derived  from  the  electrostatics  principle,  see  [15]. 

Based  on  the  polynomial  space  Pk(E),  the  discrete  discontinuous  function  spaces 
for  the  scalar  fields  and  tangential  vector  fields  are  defined  by 

=  {$  G  L°°(S)  |  VLgT:  G  Pk(E)}, 

Vjj  =  {U  G  L°°(S,  M3)  |  VLgT:  U\e  =  Ulbx  +  U2b2  with  U\  U 2  G  Pk(E )}, 

V*  =  {Ue  L°°(S,  M3)  |  VEeT  :  U\E  =  U ib1  +  U2b 2  with  Uh  U2  G  Pk(E)}. 

Because  $  G  V$,  U  G  Vu,  V  G  VJ  are  polynomials  on  each  grid  element  E,  the 
condition  <h,  U,  V  G  L°°  does  not  constitute  an  additional  constraint  to  the  discrete 
functions.  The  function  spaces  VE  and  Vj}  contain  tangential  vector  fields,  that  is 
U(x)  G  TX(S )  for  U  G  VjjUVj}.  For  every  momentum  U  G  Vu  the  restriction  U\E  is  a 
vector  field  having  polynomial  components  with  respect  to  b\  and  b2.  Test  functions 
for  the  momentum  equation  are  to  be  vector  fields  U  G  where  their  restriction 
U\E  has  polynomial  components  with  respect  to  the  dual  basis  b 1  and  b2. 

Remark  1  The  discrete  function  spaces  Vjj  and  Vfj  for  the  tangential  vector  fields 
incorporating  spherical  triangular  coordinates  ensure  the  two-dimensional  represen¬ 
tation  of  the  momentum  in  (1).  This  denotes  a  reduction  compared  to  the  three- 
dimensional  representation  in  the  Cartesian  coordinate  system  used  in  [11]  and  [13]. 
Further,  this  approach  avoids  any  kind  of  projection  step  incorporating  a  discrete 
version  of  the  Lagrangian  multiplier  in  the  numerical  scheme. 
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3.3  Space  discrete  formulation 

The  starting  point  for  the  space  discrete  formulation  is  an  approriate  integral  form 
of  the  conservation  law.  This  is  obtained  multiplying  (1)  with  a  smooth  (continuous 
in  S  with  derivatives)  test  function  p  =  (<p,  U)T ,  assuming  a  smooth  solution  q  of 
(1),  integrating  over  E  G  T  and  applying  (10),  that  is 

I  (p  ■  dtq  -  f(q )  :  Vsp)  dx  +  J  p-  f(q )  •  vEda  =  j  p-  F(x,  q)  dx. 

E  dE  E 

Here  f(q )  :  V sp  =  f®(q)  ■  V s<p  +  Y^=  1  e*  '  fu{q )  •  Vs(U  •  e*)  and  is  the  normal 
unit  vector  outward  on  dE.  This  integral  form  of  (1)  is  to  be  the  condition  that  the 
space-discrete  solution  qh(t )  G  x  VE  has  to  fulfill,  that  is 

Vp  G  14  x  V*,  VEeT 

J  (p  ■  dtqh  -  f(qh )  :  V  sp)  dx  +  J  pvn  ■  hE(x,  q™,  q°hut)  da  =  J  p  ■  F(x,  qh)  dx.  (4) 

E  dE  E 

By  means  of  the  discrete  equation,  qh  as  well  as  the  test  functions  p  are  in  discrete 
function  spaces.  The  function  space  V{j  is  used  for  V  instead  of  VE  to  simplify 
the  discrete  representation  (6)  of  (4).  Further,  due  to  the  discontinuities  of  qh  = 
(<£>,  U)  =  (<f>,  <£>tt)  along  the  edges  of  the  triangles,  the  values  of  the  flux  function 
f(qh)  are  not  defined  on  the  boundaries  dE.  That  is  why,  in  the  boundary  integral 
of  (4)  the  flux  f(qh)  is  replaced  by  the  Rusanov  numerical  flux 

=  )[(/(«n + /(O)  ■  m*)  -  %r  -  €% 

with  the  maximum  wave  speed  A  =  rnax(|iim  •  vE |  +  V 4>m,  \u°ut  ■  vE\  +  \/ <f>ottt)  in 
system  (1).  Below  the  notation  hE  =  (/i$,  hE)T  regarding  the  scalar  and  momentum 
components  will  be  used. 

To  obtain  a  matrix  formulation  of  (4),  the  decomposition  qh  =  (<£>,  U)  with  $  £  R, 
and  U  £  VE  is  regarded.  Using  the  decomposition  (3)  with  regards  to  the  Lagrange 
basis  in  E  E  T  yields 

Nk  Nk  2 

$*(*)¥>* 0*0,  U(x,t)  =  22^2uj(t)<pi(x)bi(x), 

i= 1  i=  1  1=1 

the  representation  of  qh  in  E  in  terms  of  its  component  vector 

Qh,E  =  (4*1,  »,  ®Nk,Ul,..,  Ulfk,U^,..,  U%fk),  qh  =  ( qh,E)EeT ■ 

As  expected,  the  tangential  momentum  field  U  in  E  is  represented  by  the  last  2 TV*, 
components  of  qh.E  only.  Then,  qh  is  a  solution  of  (4)  if  and  only  if  in  every  E  G  T 
for  all  j  =  1, ..,  Nf,.  and  l  —  1,  2  the  following  equations  hold 


Nk  n  p 

^  dt$i  /  Pjifidx-  /  U(qh) 

■  Vspj  dx  +  j 

J 

f  da  =  0, 

(5) 

1  J  J 

E  E 

dl 

Nk  j.  /. 

J2dtUi  /  <Pj<Pidx—  /  flu 

■  V sVj  dx  +  J 

J 

f  ipjhlLTda  =  f  (pjFfjdx. 

J 

(6) 

7-1  •'  ^ 

*  E  E 

dl 

5  E 
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Here  the  notation  f\j  =  bl  ■  fv(qh )  hlv  =  bl  ■  hv,  Flv  =  bl  ■  Fv  -  Y%m=i  /iTrL> 
F{j  =  bl  ■  Fv,  fff1  =  If  ■  fu{qh)  ■  bm  has  been  used.  Using  the  dual  momentum  space 
V{j  for  the  test  functions  in  (4)  leads  to  the  desirable  separation  of  equations  for  the 
momentum  components  for  and  (Uf)i=\,..,Nk  in  (6).  Eq.  (6),  for  l  —  1,2, 

are  two  space  discrete  momentum  equations,  only. 

To  proceed  with  the  space  discrete  system  further  below  Eqs.  (5)  and  (6)  are 
written  in  the  compact  form 


dqh 

dt 


L(qh) 


(7) 


with  an  appropriate  right  hand  side  operator  L.  The  evaluation  of  L(qfl)  includes 
the  evaluation  of  integrals  over  E  and  d E  using  the  representation  formulae  (11). 
The  integrands  include  as  well  the  flux  function  /  in  (1)  as  the  surface  geometry 
represented  by  Gram’s  determinant  g.  In  each  triangle  E  quadrature  rules  of  order 
2k  are  applied  given  in  [24],  [7],  [19],  [6].  On  each  edge  of  E  standard  Gauss-Lobatto 
rules  of  order  2k  —  1  are  applied,  since  Fekete  points  are  in  fact  Gauss-Lobatto 
points  along  the  edges,  see  [2],  Although  [5]  indicates,  that  quadrature  rules  of 
order  2k  +  1  along  the  edges  are  to  be  used  for  k  +  1-order  formal  accuracy,  this 
would  slow  down  the  scheme.  Experiments  with  a  strong  DG  method  in  [11]  have 
shown  satisfying  results  applying  2k  —  1  order  quadrature.  Furthermore,  we  rely  on 
the  superconvergence  property  of  outflow  flux  integrals  analyzed  in  [16]  to  motivate 
our  use  of  2k  —  1  quadrature  for  the  boundary  integrals. 


Remark  2  h,E{x,q™ is  a  function  of  space  and  the  conserved  variables  only, 
that  is  hE  is  independent  on  the  local  coordinate  mapping  7#.  On  the  other  hand, 
the  formulations  (5)  and  (6)  depend  on  the  coordinate  mapping  7 e  because  these 
determine  the  components  <3>j  and  U\  regarding  the  polynomial  basis  ((pi)i,..,Nk  of 
Pk(E). 


Remark  3  If  the  triangulation  T  is  time  independent,  the  mass  matrix  entries 
Mij  =  fE  (pjPi  dx  can  be  pre-evaluated  once.  Thus,  multiplying  (5)  and  (6)  by  M~l 
leads  to  the  substitution  of  ipj  and  V sTj  by  and  1 V spj ,  respectively.  The 

resulting  equations  allow  the  evaluation  of  dtqh  avoiding  any  runtime  mass  matrix 
inversion. 


3.4  Runge-Kutta  Method 

A  strong  stability-preserving  (SSP)  explicit  third-order  Runge-Kutta  (RK)  method, 
see  [14],  is  used  to  solve  the  ordinary  differential  equation  (7),  that  is  for  every  time 
step  tn  — >  tn+l 


where  s 


i—  1 


qO)  =  q(t)  =  q(0)  +  r  ^  CijL(q(j)),  i  =  1, ...,  s 


j= 0 


Qh+1  =  Q{s\ 


3, 


.  1  1  1  1  2 

(bid) C20; 0-2 1 , C30, C31, C32)  ( 


(8) 
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(9) 


Ax  [km]  /  k 

2 

4 

6 

8 

2058 

174/17920 

137/28800 

1041 

188/15360 

119/38400 

87/71680 

69/115200 

522 

94/61440 

60/153600 

44/286720 

261 

47/245760 

30/614400 

Table  1:  Depending  on  the  grid  resolution  Ax  [km]  and  the  polynomial  order  k,  the 
table  contains  (  model  resolution  h  [km]  /  number  of  grid  unknowns  ). 

SSP  methods  combine  high-order  accuracy  with  stability  properties  respecting  a 
CFL-condition.  For  nonlinear  scalar  conservation  laws  in  one  space- dimension  [22] 
showed,  that  if  the  forward  Euler  method  is  total  variation  diminishing  (TVD),  then 
(8)  is  TVD  with  an  appropriate  CFL-condition,  too.  Although  this  does  not  prove 
stability  for  (8)  applied  to  (7),  it  gives  a  good  indication,  that  this  SSP  method  is 
not  a  source  of  spurious  oscillations  in  the  discrete  solution. 

For  a  given  polynomial  degree  k  in  Sec.  3.2  a  SSP  method  of  order  k  +  1  would 
be  desirable.  Only  a  third-order  method  (9)  has  been  chosen,  because  all  SSP-RK 
methods  of  higher  order  suffer  from  more  restrictive  CFL-conditions  or  require  the 
construction  of  an  adjoint  operator.  All  numerical  experiments  in  Sec.  4  show  stable 
results  even  for  k  >  2. 


4  Numerical  Results 

The  RK-DG  method  described  in  Sec.  3  has  been  used  to  implement  a  barotropic 
model  of  the  atmosphere  which  has  been  validated  performing  numerical  experi¬ 
ments.  The  validation  process  is  carried  out  in  two  steps.  At  first,  a  convergence 
study  considering  steady-state  and  unsteady  analytical  solutions  of  the  nonlinear 
SWE  is  performed.  After  that  a  barotropic  instability  in  a  localized  jet  stream  has 
been  carried  out. 

For  the  test  with  an  analytical  solution  q  =  (<f>ana,  ^anaUana)T ^  the  numerical  error 
is  evaluated  using  the  normalized  L2-error  of  the  geopotential  held  <P/)  ,  that  is 

UK  \  ll'hana  —  <h/i || L2(S) 

"(#k)  =  ' 

For  every  element  E  a  local  grid  resolution  A xe,  the  size  of  the  largest  edge  in 
E ,  and  a  local  model  resolution  He  is  defined  by  He  =  y/\E\/Nk-  Grid  resolution 
Ax  and  model  resolution  h  arise  from  the  corresponding  maximum  values  over  all 
elements.  Table  1  contains  h  and  the  number  of  grid  unknowns,  depending  on  Ax 
and  the  polynomial  order  k. 

As  a  result  of  numerical  experiments  in  every  element  E  the  CFL-condition 

AtE  3 

XE- — =  A(1  lV  At  <  AtcFL  =  mm AtE 
hE  4(fc  +  1)  egt 

was  chosen  where  A^  is  the  characteristic  wave  speed  A^  =  |u|  +  \/$.  All  numerical 
experiments  gave  stable  results  without  spurious  oscillations  and  a  relaxed  condition 
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Figure  2:  Section  4.1  (Steady-state  Solid  Body  Rotation),  Normalized  L2-error  rj($h) 
after  integration  time  5  days  as  a  function  of  model  resolution  h,  k  =  2  (■■•),  k  —  4 
(— ),*  =  6(-),*  =  8(-). 


showed  numerical  instabilities.  In  the  case  k  =  2,  with  third  order  accuracy  in  space 
and  time,  this  condition  simplifies  to 


A tE  _  1 
hE  4’ 


which  is  similar  to  the  CFL-condition  for  linear  stability  in  [5]  in  the  one-dimensional 
scalar  case.  The  critical  time  step  length  A tcFL  is  derived  in  every  time  step,  for 
given  values  hE  and  \E.  All  numerical  experiments,  except  Sec.  4.2,  have  been 
performed  with  the  time  step  At  =  AtcFL-  Due  to  accuracy  limitations  the  smaller 
time  step  At  =  At(^FL  is  used  for  the  experiments  in  Sec.  4.2.  With  this  time 
step  control  the  model  shows  stable  numerical  results  and  uses  reasonable  time  step 
length. 


4.1  Steady-state  Solid  Body  Rotation 

This  test  contains  a  steady-state  solution  of  the  nonlinear  SWE,  see  [26,  case  2],  The 
velocity  held  u  is  a  westerly  wind  with  the  meridional  distribution  of  a  solid  body 
rotation.  The  geopotential  height  $  is  given  in  geostrophic  balance  to  u.  Thus,  for 
the  duration  of  the  integration  the  initial  data  have  to  be  maintained. 

For  the  validation  [26]  recommend  the  evaluation  of  the  normalized  L2-error  //($/,  ) 
after  an  integration  time  of  5  days.  Fig.  2  shows  rj(<&h)  for  different  polynomial 
orders  k  =  2,4,  6,8  as  a  function  of  the  model  resolution  h.  For  all  choices  of  the 
polynomial  order  k  the  model  converges  and  reduces  the  error  almost  up  to  machine 
precision  for  k  =  6,  8.  The  errors  decrease  significantly  for  increasing  k.  Table  2 
shows  the  expected  order  of  convergence  k  +  1.  These  results  are  very  close  to  the 
convergence  studies  in  [12],  [20]  and  [11].  For  this  steady-state  solution  no  limitation 
due  to  the  third  order  RK  method  is  observable. 
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Test  /  k 

2 

4 

6 

8 

Section  4.1 

2.86 

4.97 

6.97 

8.78 

Section  4.2 

3.02 

5.02 

4.83 

4.20 

Section  4.3 

4.84 

5.70 

7.25 

8.79 

Table  2:  Sections  4.1,  4.2,  4.3,  Experimental  order  of  convergence  for  k  =  2,4,  6,8. 


Figure  3:  Section  4.2  (Unsteady  Solid  Body  Rotation),  Normalized  L2-error  v(^h) 
after  integration  time  5  days  as  a  function  of  model  resolution  h,  k  =  2  k  =  4 

(— ),*  =  6(-),*  =  8(-). 

4.2  Unsteady  Solid  Body  Rotation 

This  test  concerns  an  unsteady  solution  of  the  nonlinear  SWE,  see  [17,  example  3]. 
Similar  to  the  last  test,  the  velocity  field  u  is  a  solid  body  rotation,  but  with  an 
inclination  to  the  Earth’s  rotation  axis.  An  adequate  unsteady  geopotential  field 
is  available  and  an  axially  symmetric  orographic  field  has  to  be  regarded.  This 
analytical  solution  moves  into  a  westerly  direction  and  has  a  time  period  of  one  day. 

Again,  the  normalized  L2-error  is  evaluated  after  an  integration  time  of 

5  days.  In  this  Section  the  time  step  is  chosen  as  At  =  At(^FL ;  which  is  different 
from  all  other  numerical  experiments  in  this  article.  Fig.  3  shows  for  different 

polynomial  orders  k  =  2,4,  6,8  as  a  function  of  the  model  resolution  h.  For  all 
choices  of  the  polynomial  order  k  the  model  converges,  even  in  this  unsteady  test 
case.  The  errors  decrease  significantly  for  increasing  k. 

For  all  experiments  with  small  errors  (r)(&h)  <  1CU9)  the  limiting  factor  for 
accuracy  is  the  time  step  At  and  not  the  model  resolution  h  any  more.  To  see  this, 
Fig.  4  is  given,  which  shows  ^(4>/))  as  a  function  of  At,  for  fixed  parameters  ( h ,  k). 
All  experiments  which  are  not  plotted  in  Fig.  4  give  the  same  results  independent 
of  the  choice  of  At  <  AtcFL ■  It  is  to  be  seen,  that  r](&h)  stagnates  for  decreasing 
h ,  for  fixed  At  =  AtcFL  and  all  experiments  with  small  errors.  At  the  same  time, 
for  these  experiments  Fig.  4  shows  0(At3)  convergence  of  r)(&h)-  Both  observations 
together  yield,  once  r](&h)  is  small  enough,  only  a  decreasing  time  step  leads  to 
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Figure  4:  Section  4.2  (Unsteady  Solid  Body  Rotation),  Normalized  L2-error  r)($h) 
after  integration  time  5  days  as  a  function  of  time  step  At/ AtCFL,  (h,  k)  =  (60km,  4) 

(• - ),  (h,k)  =  (30km,  4)  (A - ),  (h,k)  =  (87km,  6)  (•  —  •),  (h,k)  =  (44km,  6) 

(A  —  ■),  (h,k)  =  (137km,  8)  (•—),.  (h,k)  =  (69km,  8)  (A—). 

decreasing  A  further  error  reduction  is  anticipated  for  smaller  time  steps 

than  ,  which  is  not  shown  here  due  to  limited  computational  resources.  Of 

course  this  effect  could  be  avoided  using  an  appropriate  RK  method  of  order  k  +  1. 
As  explained  in  Sec.  3.4  a  SSP-RK  method  with  a  tolerant  CFL-condition  is  not 
available. 

Table  2  shows  for  k  =  2,4  the  expected  order  of  convergence  k  +  1.  Because  for 
k  =  6,8  the  convergence  is  limited  by  the  third  order  time  step,  in  this  case  the 
convergence  rates  in  space  are  suboptimal.  At  the  same  time  the  expected  third 
order  accuracy  in  time  can  be  achieved.  Because  the  absolute  L2-error  ||<hana||z,2(S) 
is  larger  compared  to  the  similar  steady-state  case  in  Sec.  4.1,  for  the  polynomial 
orders  k  —  2,4,6  the  normalized  error  is  even  smaller.  The  error  for  k  =  8  is  limited 
by  the  time  step,  which  results  in  slightly  larger  errors  compared  to  Sec.  4.1. 

4.3  Unsteady  Jet  Stream 

This  test  contains  a  second  unsteady  solution  of  the  nonlinear  SWE,  see  [17,  example 
4],  This  time,  the  velocity  held  u  is  an  axially  symmetric  westerly  wind  jet  stream 
superimposed  by  a  smaller  solid  body  rotation,  where  the  axis  of  the  jet  stream 
is  inclined  to  the  Earth’s  rotation  axis.  Due  to  the  jet  stream  strong  meridional 
gradients  are  present,  which  presents  an  additional  difficulty  compared  to  Sec.  4.2. 
As  in  Sec.  4.2,  an  adequate  unsteady  geopotential  held  is  available  and  an  axially 
symmetric  orographic  held  has  to  be  regarded.  The  solution  moves  into  westerly 
direction  and  has  a  time  period  of  one  day. 

The  normalized  L2-error  r)(&h)  is  evaluated  after  an  integration  time  of  5  days. 
Fig.  5  shows  the  normalized  L2-error  r](<&h)  for  different  polynomial  orders  k  = 
2,4,  6,8  as  a  function  of  the  model  resolution  h.  The  method  shows  experimental 
convergence  for  this  test  case  with  strong  meridional  gradients.  For  increasing  poly- 
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Figure  5:  Section  4.3  (Unsteady  Jet  Stream),  Normalized  L2-error  after  inte¬ 
gration  time  5  days  as  a  function  of  model  resolution  h,  k  =  2  (•••),  k  =  4  ( - ), 

k  —  6  (— •),  k  —  8  (—). 

nornial  order  k  the  error  rj^h)  decreases,  but  remains  significantly  larger  than  in 
Sec.  4.2.  This  is  probably  the  reason  for  the  order  of  convergence  k+  1,  see  Table  2, 
without  any  limitation  due  to  the  third  order  RK  method  in  this  unsteady  solution. 
Due  to  the  locality  of  the  jet  stream,  the  9th-order  experiment  ( k  =  8)  on  a  coarser 
grid  resolves  this  test  not  as  well  as  the  7th-order  ( k  =  6)  experiment  which  explains 
the  crossing  of  both  error  lines. 

4.4  Perturbed  Jet  Stream 

Supplementing  the  standard  tests  for  atmospheric  models  based  on  SWE,  [9]  pro¬ 
posed  a  test  describing  a  barotropic  instability.  The  initial  velocity  held  u  is  an 
axially  symmetric  westerly  jet  stream  with  the  same  axis  as  the  Earth’s  rotation 
axis.  As  in  Sec.  4.3,  this  jet  stream  includes  strong  meridional  gradients  and  con¬ 
stitutes  a  rather  local  feature.  Based  on  u  a  geopotential  height  <f>  is  derived  in 
geostrophic  balance  to  u.  Additionally  a  small  perturbation  <f>(2  is  added,  such  that 
the  initial  condition  for  the  test  is  the  geopotential  held  $  +  JD-  As  a  consequence  of 
the  perturbation  $,2,  this  experiment  should  not  maintain  the  initial  data.  [9]  give 
a  detailed  description  of  the  barotropic  instability  developing  within  the  jet  stream 
from  day  four  to  day  six. 

At  hrst,  the  model  should  be  able  to  maintain  the  geostrophic  balanced  how,  as 
long  as  the  initial  perturbation  does  not  lead  to  instabilities.  Unlike  for  spectral 
models,  this  is  not  trivial  for  a  model  with  a  grid  which  is  not  aligned  to  the  zonal 
how.  Fig.  6  shows  the  vorticity  held  after  four  days  in  a  cutout  of  the  global  model 
domain  for  k  =  2  and  two  different  model  resolutions.  Whereas  a  pronounced  zonal 
wavenumber  hve  is  visible  for  h  =  94km  the  experiment  with  the  higher  resolution 
h  =  47km  reduces  the  grid  effect  considerably.  The  same  plot  after  four  days  is  given 
in  the  upper  plot  of  Fig.  7  for  the  model  resoltuion  h  =  87km  but  polynomial  order 
k  =  6.  This  “high  order“  experiment  seems  to  be  more  adequate  to  maintain  the 
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Latitude  Latitude 


Figure  6:  Section  4.4  (Perturbed  Jet  Stream),  Vorticity  (contour  interval  2  x 
10~5s"1),  day  4,  k  —  2,  top:  h  =  94km,  bottom:  h  =  47km. 
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Figure  7:  Section  4.4  (Perturbed  Jet  Stream),  Vorticity  (contour  interval  2  x 
10_5s_1),  k  =  6,  top:  day  4,  h  =  87km,  middle:  day  6,  h  =  87km,  bottom: 
day  6,  h  =  44km. 
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Figure  8:  Section  4.4  (Perturbed  Jet  Stream),  Vorticity  (contour  interval  2  x 
10_5s_1),  (h,k)  =  (30km,  4),  top:  day  4,  bottom:  day  6. 

geostrophic  structure  of  the  flow  avoiding  a  low  frequency  zonal  wave.  The  two  lower 
plots  of  Fig.  7  show  the  vorticity  for  k  =  6  and  the  resolutions  87km  and  44km  after 
six  days.  In  both  experiments  the  laminar  flow  structure  looks  very  similar,  which 
gives  an  indication  for  the  experimental  convergence  in  smooth  regions.  In  contrast 
to  this,  the  areas  with  strong  gradients,  where  barotropic  instabilities  develop,  show 
spurious  oscillations.  The  same  experiment  has  been  performed  with  polynomial 
order  k  =  4  and  high  model  resolution  h  =  30km,  see  Fig.  8.  Compared  to  the 
“low“  and  “high“  order  experiments  in  Figs.  6  and  7  the  geostrophic  balanced 
flow  is  maintained  as  well  as  the  barotropic  instabilities  show  only  small  spurious 
oscillations.  Fig.  8  shows  very  good  agreement  with  the  inviscid  run  of  [9]  and  the 
results  in  [23]. 

All  numerical  experiments  have  shown  the  discrete  global  conservation  of  mass 
up  to  machine  precision,  which  is  a  direct  consequence  of  the  conservation  properties 
of  the  DG  method.  The  discrete  global  conservation  of  energy  Eh  =  +  (JA  + 

J?#)2  —  &2B  cannot  be  anticipated,  because  Eh  is  not  a  conserved  variable  of  the 
hyperbolic  system  (1).  To  investigate  the  discrete  energy  conservation,  Fig.  9  shows 
the  absolute  value  of  the  normalized  energy  error 

f  T?  \  fsM  —  ^ana)  d X 

v(Eh)  =  PjAs) 
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Figure  9:  Section  4.4  (Perturbed  Jet  Stream),  normalized  energy  error  \r)(Eh)\  as  a 
function  of  model  resolution  h,  k  —  2  (•  •  •),  k  —  4  ( - ),  k  =  6  (— •),  k  =  8  (— ). 

as  a  function  of  the  model  resolution  h.  The  energy  conservation  is  better  than 
2.0  x  10-4,  even  for  the  coarse  model  resolution  of  180km.  Furthermore,  for  fixed 
polynomial  orders  k  and  decreasing  h  the  energy  error  convergences.  In  the  same 
manner  for  a  fixed  grid  resolution  Ax,  an  increasing  k  leads  to  an  energy  error 
convergence.  The  time  evolution  of  r)(Eh)  is  plotted  in  Fig.  10.  For  each  k  only  the 
experiment  with  the  smallest  model  resolution  h  is  shown,  see  Table  1.  For  almost 
all  experiments  an  energy  loss  is  observed.  For  k  =  2  the  error  is  much  larger  than 
for  k  >  4.  This  could  be  an  indication  for  a  loss  of  accuracy  due  to  the  2k  —  1-order 
quadrature  rule  along  the  element  edges.  As  long  as  the  geostrophic  flow  remains 
stable,  up  to  day  four,  the  energy  error  remains  very  small.  After  the  barotropic 
instability  develops,  after  day  four,  all  experiments  show  an  enforced  loss  of  energy. 
To  explain  this  phenomenon  we  find,  that  the  flow  in  the  instability  regions  tends 
to  develop  discontinuities,  which  results  in  bigger  jumps  in  the  discrete  solution. 
Considering  this,  [4]  describes  for  a  linear  example,  that  the  rate  of  dissipation  in  the 
DG  method  is  given  by  the  jumps  of  the  solution.  This  property  gives  an  explanation 
for  the  stronger  energy  loss  during  the  development  of  barotropic  instabilities. 

5  Summary  and  Outlook 

A  global  barotropic  model  of  the  atmosphere  has  been  developed  based  on  the 
spherical  shallow  water  equations  and  discretized  by  a  Runge-Kutta  discontinuous 
Galerkin  method.  The  equations  are  formulated  as  a  hyperbolic  conservation  law 
on  the  sphere,  a  two-dimensional  surface  in  M3,  using  coordinate  independent  dif¬ 
ferential  opertators. 

Spherical  triangular  coordinates,  which  are  local  coordinate  mappings  7 e,  and 
high-order  polynomial  spaces  are  defined  on  each  curved  element  E  C  S  of  a  given 
spherical  triangluar  grid.  7#  yields  a  two-dimensional  representation  of  the  tangen¬ 
tial  momentum  fields  and  only  two  momentum  equations  which  denotes  a  reduction 
of  variables  compared  to  the  three-dimensional  representation  in  a  Cartesian  coordi- 
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Figure  10:  Section  4.4  (Perturbed  Jet  Stream),  normalized  energy  error  r)(Eh)  as  a 

function  of  integration  time  T,  ( h,k )  =  (47km,  2)  (•  •  •),  (h,k)  =  (30km,  4)  ( - ), 

(h,k)  =  (44km,  6)  ( — ),  (h,k)  =  (69km,  8)  (— ). 

nate  system.  Governed  by  an  integral  form  of  the  conservation  law  the  space  discrete 
discontinuous  Galerkin  method  is  obtained  including  a  Rusanov  numerical  flux.  An 
explicit  third  order  Runge-Kutta  method  (strong  stability-preserving)  is  applied  to 
the  space  discrete  system.  The  discrete  method  avoids  any  kind  of  explicit  smooth¬ 
ing  such  as  diffusion  or  filter  operators.  A  barotropic  model  of  the  atmosphere  is 
implemented  based  on  the  Runge-Kutta  discontinuous  Galerkin  method  completed 
by  a  time  step  control  based  on  a  CFL-condition.  The  model  is  validated  regarding 
numerical  experiments. 

Steady-state  and  unsteady  analytical  solutions  of  the  nonlinear  shallow  water 
equations  are  considered.  In  all  cases  the  model  achieves  experimental  convergence 
to  the  known  solution.  Furthermore,  for  increasing  polynomial  order  k,  the  error 
decreases  significantly.  For  the  steady-state  solid  body  rotation  and  the  unsteady 
jet  stream  the  experimental  orders  of  convergence  are  k  +  1.  This  convergence  rate 
is  achieved  for  the  unsteady  solid  body  rotation,  only  if  k  <  4.  For  higher  orders 
k  and  very  small  errors  the  time  step  At  limits  the  error  (third  order  Runge-Kutta 
method).  This  result  of  the  presented  atmospheric  model  demonstrates,  that  the 
model  time  step  is  restricted  not  only  by  a  CFL-condition  but  by  accuracy  demands, 
too.  Further,  a  qualitative  test  regarding  a  barotropic  instability  within  a  perturbed 
jet  stream  is  performed.  In  very  good  agreement  to  published  results  the  experiment 
shows  the  evolution  of  a  breaking  barotropic  wave  after  six  days  without  generating 
spurious  oscillations.  The  global  energy  error  converges  with  decreasing  grid  resolu¬ 
tion  as  well  as  with  increasing  polynomial  order.  For  all  numerical  experiments  the 
anticipated  mass  conservation  is  achieved  up  to  machine  precision. 

The  presented  results  show  the  great  potential  of  the  Runge-Kutta  discontinuous 
Galerkin  method  in  a  simplified  atmospheric  model.  Thus,  the  application  to  three 
dimensional  atmospheric  equations,  which  form  a  hyperbolic  system,  seems  to  be 
warranted.  Further,  the  method  demonstrates  the  application  of  local  coordinate 
mappings,  here  spherical  triangular  coordinates,  and  its  capacity  to  represent  tan- 
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gential  momentum  fields  that  are  locally  two-dimensional.  Using  the  same  technique 
for  global  atmospheric  models  in  three-dimensional  prismatic  grids  yields  the  useful 
decomposition  into  tangential  and  vertical  components. 
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Appendix,  Manifolds 


The  spherical  domain  S  is  regarded  as  a  two-dimensional  C'1-manifold  embedded 
in  M3.  Thus,  the  operators  and  Gaufi’  theorem  on  S  are  special  cases  of  the  stan¬ 
dard  definitions  from  differential  geometry.  Because  the  local  representation  of  the 
operators  and  the  integrals  regarding  a  coordinate  mapping  7  are  essential  to  the 
implementation  of  the  model,  these  formulas  are  given  below. 

Let  /  G  be  a  scalar  function  and  u  G  G1(S',  M3)  a  vector  field  with  u(x)  G 

TX(S).  Further,  for  a  fixed  x  G  S  let  (77,  r2)  be  a  orthonormal  basis  of  TX(S).  Then 
the  spherical  gradient  and  spherical  divergence  are  defined  by 

2  2 

V5/(x)  =  ^  d T.f(x)  Ti ,  div5  u(x)  =  Ti  ■  dTiu(x). 

i= 1  i= 1 

Then,  for  an  open  subset  E  C  S  with  smooth  boundary  a  formulation  of  Gaufi’ 
theorem  on  S  is 


j  f  div,s  u  dx 


j  u  ■  \/sfdx  + 


fu  ■  vE  da. 


E 


E 


dE 


(10) 


Now,  let  O  C  S  be  open  with  a  local  coordinate  mapping  7  :  U  C  M2  — »  U, 
D  C  U,  E  C  U  and  7 (D)  =  E.  Further,  remember  ( b 1,b2)  the  dual  basis  of  TX(S ) 
and  Gram’s  determinant  g. 

Let  D  C  M2  open,  E  C  S  open  in  S  and  7:h->£a  local  coordinate  mapping. 
7  generates  the  local  coordinate  system 


bi  =  dyi  7,  2  =  1,2 

which  is  a  basis  of  the  tangential  space  TX(S).  The  condition  bt  ■  &  =  Sj  deter¬ 
mines  the  dual  basis  (fe1^2).  Gram’s  matrix  (metric  tensor),  its  inverse  and  Gram’s 
determinant  are  defined  by 

9ij  =  bt  ■  bj ,  glJ  =  bl  ■  bJ ,  g  =  det(gij)ij=i!2- 

Derivatives  of  the  basis  functions  bi  are  expressed  using  Christoffel  symbols  which 
fulfill  the  equations 

nnk 

dy3  //  =  -Y)kbk ,  r %  =  —  (dyi gjn  +  dyj gin  -  dVn gt] ) . 
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Then,  for  y  G  D  the  operators  fulfill  the  local  representations 


Vs/ |7(y)  =  YlfivM  °'y^vb'’  div^Ml7(2/)  =  ^2  ^-[dyiVaiu  O  -b%. 

i= 1  V  ■ 9 


i=  1 


The  face  and  line  integrals  satisfy  the  local  representations 


J  f  dx  =  J  yjgf  o^fdy,  J  f  da  =  J  y/g\(b  ,  b  )  ■  vD\f  o  7  da  (11) 

ED  dE  dD 

where  uo  is  the  unit  normal  on  dD. 
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